{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VMLS\n",
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 5 \n",
    "\n",
    "# Linear independence \n",
    "\n",
    "\n",
    "### 5.1 Linear dependence \n",
    "\n",
    "### 5.2 Basis \n",
    "\n",
    "**Cash flow replication.** Let’s consider cash flows over 3 periods, given by 3-vectors. \n",
    "We know from VMLS page [93](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section*.108) that the vectors \n",
    "\n",
    "$$\n",
    "e_1 = \\begin{bmatrix} 1\\\\ 0\\\\ 1\\end{bmatrix}, \n",
    "l_1 = \\begin{bmatrix} 1\\\\ -(1+r)\\\\ 0 \\end{bmatrix},\n",
    "l_2 = \\begin{bmatrix} 0\\\\ 1\\\\ −(1 + r)  \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "form a basis, where $r$ is the (positive) per-period interest rate. The first vector $e_1$ \n",
    "is a single payment of $\\$1$ in period (time) t = 1. The second vector $l_1$ is loan of $\\$1$ \n",
    "in period $t = 1$, paid back in period $t = 2$ with interest $r$. The third vector $l_2$ is \n",
    "loan of $\\$1$ in period $t = 2$, paid back in period $t = 3$ with interest $r$. Let’s use this \n",
    "basis to replicate the cash flow $c = (1, 2,−3)$ as \n",
    "\n",
    "$$\n",
    "c = α_1e_1 + α_2l_1 + α_3l_2 = α_1 \\begin{bmatrix} 1 \\\\ 0 \\\\ 0\\end{bmatrix} + α_2 \\begin{bmatrix} 1\\\\−(1 + r) \\\\ 0 \\end{bmatrix} + α_3 \\begin{bmatrix} 0 \\\\ 1\\\\  −(1 + r) \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "From the third component we have \n",
    "$$c_3 = α_3(−(1 + r)), so$$ \n",
    "$$α_3 = \\frac{−c_3}{(1 + r)}.$$ \n",
    "From the second component we have \n",
    "$$c_2 = α_2(−(1 + r)) + α_3 = α_2(−(1 + r)) − \\frac{c_3}{(1 + r)}, so$$\n",
    "$$α_2 = \\frac{−c_2}{(1 + r)} − \\frac{c_3}{(1 + r)^2}.$$\n",
    "Finally from $c_1 = α_1 + α_2$, we have\n",
    "$$α_1 = c_1 + \\frac{c_2}{(1 + r)} + \\frac{c_3}{(1 + r)^2},$$ which is the net present value (NPV) of the cash flow $c$. \n",
    "\n",
    "Let’s check this in Julia using an interest rate of $5\\%$ per period, and the specific cash flow `c = (1, 2,−3)`. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Tuple{String,Float64},1}:\n",
       " (\"a1:\", 0.18367346938775508)\n",
       " (\"a2:\", 0.8163265306122449) \n",
       " (\"a3:\", 2.857142857142857)  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = 0.05\n",
    "e1 = [1,0,0]; l1 = [1,-(1+r),0]; l2 = [0,1,-(1+r)];\n",
    "c = [1,2,-3];\n",
    "# Coefficients of expansion\n",
    "alpha3 = -c[3]/(1+r);\n",
    "alpha2 = -c[2]/(1+r) -c[3]/(1+r)^2;\n",
    "alpha1 = c[1] + c[2]/(1+r) + c[3]/(1+r)^2 # NPV of cash flow\n",
    "\n",
    "[(\"a1:\", alpha1),\n",
    "(\"a2:\", alpha2), \n",
    "(\"a3:\", alpha3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  1.0\n",
       "  2.0\n",
       " -3.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha1*e1 + alpha2*l1 + alpha3*l2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Later in the course we’ll use an automated and simple way to find the coefficients in the expansion of a vector in a basis.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3 Orthonormal vectors\n",
    "**Expansion in an orthonormal basis.** Let’s check that the vectors \n",
    "$$\n",
    "a_1 = \\begin{bmatrix} 0\\\\ 0\\\\ -1\\end{bmatrix}, \n",
    "a_2 = \\frac{1}{\\sqrt{2}}\\begin{bmatrix} 1\\\\ 1\\\\ 0 \\end{bmatrix},\n",
    "a_3 = \\frac{1}{\\sqrt{2}}\\begin{bmatrix} 1\\\\ -1\\\\ 0\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "form an orthonormal basis, and check the expansion of $x = (1, 2, 3)$ in this basis,\n",
    "\n",
    "$$x = (a^{T}_{1} x)a_1 + · · ·+ (a^{T}_{n}x)a_n.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.0, 0.9999999999999999, 0.9999999999999999)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a1 = [0,0,-1]; a2 = [1,1,0]/sqrt(2); a3 = [1,-1,0]/sqrt(2);\n",
    "norm(a1), norm(a2), norm(a3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 0.0, 0.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a1'*a2, a1'*a3, a2'*a3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Int64,1}:\n",
       " 1\n",
       " 2\n",
       " 3"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [1,2,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.9999999999999999\n",
       " 1.9999999999999996\n",
       " 3.0               "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get coefficients of x in orthonormal basis\n",
    "beta1 = a1'*x; beta2 = a2'*x; beta3 = a3'*x\n",
    "# Expansion of x in basis\n",
    "xexp = beta1*a1 + beta2*a2 + beta3*a3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4 Gram–Schmidt algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is a Julia implementation of Algorithm [5.1](https://web.stanford.edu/~boyd/vmls/vmls.pdf#algorithmctr.5.1) in VMLS (Gram–Schmidt algorithm). It takes as input an array `[ a[1], a[2], ..., a[k] ]`, containing the $k$ vectors $a_1, . . . , a_k$. If the vectors are linearly independent, it returns an array `[ q[1], ..., q[k] ]` with the orthonormal set of vectors computed by the Gram – Schmidt algorithm. If the vectors are linearly dependent and the Gram–Schnidt algorithm terminates early in iteration `i`, it returns the array `[ q[1], ..., q[i] ]` of length `i`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gram_schmidt (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function gram_schmidt(a; tol = 1e-10)\n",
    "    \n",
    "    q = []\n",
    "    for i = 1:length(a)\n",
    "        qtilde = a[i]\n",
    "        for j = 1:i-1\n",
    "            qtilde -= (q[j]'*a[i]) * q[j]\n",
    "        end\n",
    "        if norm(qtilde) < tol\n",
    "            println(\"Vectors are linearly dependent.\")\n",
    "            return q\n",
    "        end\n",
    "        push!(q, qtilde/norm(qtilde))\n",
    "        end;\n",
    "    return q\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On `line 3`, we initialize the output array as the empty array. In each iteration, we add the next vector to the array using the push! function (`line 13`).\n",
    "\n",
    "**Example.** We apply the function to the example on page [100](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section*.117) of VMLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Array{Int64,1},1}:\n",
       " [-1, 1, -1, 1]\n",
       " [-1, 3, -1, 3]\n",
       " [1, 3, 5, 7]  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [ [-1, 1, -1, 1], [-1, 3, -1, 3], [1, 3, 5, 7] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Any,1}:\n",
       " [-0.5, 0.5, -0.5, 0.5]\n",
       " [0.5, 0.5, 0.5, 0.5]  \n",
       " [-0.5, -0.5, 0.5, 0.5]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = gram_schmidt(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Array{Tuple{String,Float64},1}:\n",
       " (\"norm(q[1]:)\", 1.0)\n",
       " (\"q[1]'*q[2]:\", 0.0)\n",
       " (\"q[1]'*q[3]:\", 0.0)\n",
       " (\"norm(q[2]):\", 1.0)\n",
       " (\"q[2]'*q[3]:\", 0.0)\n",
       " (\"norm(q[3]):\", 1.0)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# test orthnormality\n",
    "[(\"norm(q[1]:)\",norm(q[1])),\n",
    "(\"q[1]'*q[2]:\",q[1]'*q[2]),\n",
    "(\"q[1]'*q[3]:\",q[1]'*q[3]),\n",
    "(\"norm(q[2]):\",norm(q[2])),\n",
    "(\"q[2]'*q[3]:\",q[2]'*q[3]),\n",
    "(\"norm(q[3]):\",norm(q[3]))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example of early termination.** If we replace $a_3$ with a linear combination of $a_1$, and $a_2$ the set becomes linearly dependent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Array{Float64,1},1}:\n",
       " [-1.0, 1.0, -1.0, 1.0]\n",
       " [-1.0, 3.0, -1.0, 3.0]\n",
       " [-1.8, 2.8, -1.8, 2.8]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = [ a[1], a[2], 1.3*a[1] + 0.5*a[2] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vectors are linearly dependent.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2-element Array{Any,1}:\n",
       " [-0.5, 0.5, -0.5, 0.5]\n",
       " [0.5, 0.5, 0.5, 0.5]  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = gram_schmidt(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example of independence-dimension inequality.** We know that any three $2$ - vectors must be dependent. Let’s use the Gram-Schmidt algorithm to verify this for three specific vectors.',"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vectors are linearly dependent.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2-element Array{Any,1}:\n",
       " [0.707107, 0.707107] \n",
       " [-0.707107, 0.707107]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "three_two_vectors = [ [1,1], [1,2], [-1,1] ];\n",
    "q = gram_schmidt(three_two_vectors)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.1",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
